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Abstract 

A classic problem, the design of the tallest column, is solved again using a different 
method. By the use of a similarity solution the equations are transformed and the 
difficult singularity at the cndpoint is peeled away. The resulting autonomous system 
has a critical point and the solution must be on its stable manifold. The solution 
is found by starting near the critical point in the direction of the stable manifold, 
and solving backwards numerically This removes the need for an iterative integration 
method that was previously used. The method is shown to work for clamped or hinged 
boundary condition and can also be used for other problems involving singularities at 
the endpoints. 

Keywords: Tallest column, eigenvalue optimization, singular ODE, similarity solution, sta- 
ble manifold. 



1 Introduction 

We reexamine a classic eigenvalue maximization problem: The design of the tallest column. 
This problem has a lineage that originates from Euler's analysis of column buckling. In 
[1] Keller determines the tapering of a thin column with fixed volume which maximizes its 
buckling load. Later, in [5], Keller and Niordson solved the problem of designing the tallest 
free standing column. The equations have a nasty singularity at the top of the column. In 
[5] this singularity is neutralized by formulating an equivalent integral equation. The latter 
is solved by numerical iteration. 

More recently, in [3] and [7j, Cox and McCarthy showed that the operator involved can 
have a continuous spectrum and they find the tapering of the tallest column using other 
methods. In this paper, no claims are made as to the validity of the equations involved 
or the derivation of them, see [6], [8], [3] and [7] for a discussion regarding the existence 



*Department of Mathematics, UC Berkeley, Berkeley CA; e-mail: yfarjoun@math.berkeley.edu 
^Department of Mathematics, UC Berkeley, Berkeley CA 



1 



of the optimal design and the control requirements of the design. Nevertheless, Keller and 
Niordson's solution is highly plausible since the problem is infinitesimally close to a problem 
with a discrete spectrum. 

The current paper concentrates on solving the difficult equations ()19l — 1230 . originally derived 
by Keller and Niordson, by a new method. 

First, in Section^ we derive the boundary value problem (BVP) which is an ordinary differ- 
ential equation (ODE) system for the deflection of the column and the cross-sectional area, 
together with boundary conditions (BC). The ODE's contains the buckling load eigenvalue, 
to be determined as part of the solution. This derivation follows the derivation in [5] closely 
and is included for completeness. 

In Section [3] we solve the BVP. First we show that the ODE's have a scaling symmetry 
and an exact similarity solution which has the desired asymptotic behavior at the tip of the 
column. This is where the singularity is. The similarity solution suggests a transformation 
of variables that "peels away the singularity" at the tip. This transformation yields ODE's 
for the new variables which form an autonomous system (AS). In the new variables, the 
similarity solution is represented by a critical point of the AS. Also, due to another symme- 
try, the eigenvalue disappears from the equations and now appears only in the BC at the 
base of the column. The transformed BC dictate that the solution must start on a certain 
surface in the phase space of the AS and approach the critical point. 

Since the critical point of the AS has both stable and unstable manifolds, solving the AS 
numerically from the base to the critical point is difficult — the solution tends to deviate 
to a growing solution. Instead, we examine the behavior of the AS near the critical point 
in order to solve the system backwards, from the tip (the critical point) to the base. By 
linearizing the AS around the critical point, we identify stable and unstable directions. 
Starting near the critical point, on the stable manifold, the AS is solved numerically, until 
the solution intersects the surface that the BC defines. The point of intersection with the 
surface determines the eigenvalue and the problem is solved. The numerical results agree 
with those of Keller and Niordson. 

In Section m the problem is solved for other boundary conditions. This is inexpensive once 
all the groundwork has been done. 

This new solution method applies to other eigenvalue maximization problems whose equa- 
tions have certain scale-invariant structure. In a companion paper this method determines 
the tapering of a javelin so the its first mode has the highest frequency of vibration (subject 
to length and volume constraint). 
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2 Derivation of the Boundary Value Problem 



2.1 Setup 

Consider columns, all of the same volume, clamped at the base and free at the top. We 
want to find the shape of the tallest column that will not buckle under its own weight. As 
in the Keller and Niordson papers, we solve the problem for a specific class of permissible 
designs. We assume that the column is thin, i.e. the characteristic width is much less than 
the height of the column. In addition, we only allow columns with geometrically similar, 
equally oriented and convex cross-sections. See figure ([1]). 




Figure 1: The shape of the column is governed by the cross-sectional area function a(s) 
measuring the area of a cross-section at a point located at arclength s measured along the 
center axis of the column from its tip. All the cross-sections are convex, are geometrically 
similar, and are equally oriented. 

We parameterize the column by arclength s, measured from the tip of the column, along its 
center axis. The design information is contained in a single function a{s), the cross-sectional 
area at s. Lastly, we concern ourselves with the bending of the column in a specified plane 
only. This allows us to specify the configuration of the column using a single function 
9{s) — the angle that the center axis of the column makes with the vertical, measured at 
point s. See Figure ([2]). 
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z(s) 



Figure 2: The shape of the bent column is specified by the angle that the center axis of the 
column makes with the vertical at each point. The points on the column are parameterized 
by arclength s, measured from the tip along the columns center axis. z{s) is the vertical 
height of the point s. 



To find the BVP satisfied by 9[s) for a given cross-section a(s), we write the total energy 
of the system as a functional of and compute its variational equations. 

2.2 Energy Minimization 

The total energy of the column, elastic and gravitational, is a functional of 9: 



Here, L is the length of the column, g is the gravitational acceleration, p is the density of the 
material, and z{s) is the vertical elevation at s. The bending modulus, h{s) is proportional 
to a^(s), 



Here, c is a dimensionless constant determined by the cross-sectional shape and orientation, 
and E is Young's modulus. For a discussion of the choice of cross-sectional designs and 
derivation of the bending modulus, see [4]. Since z{s) is the vertical elevation of a particle 




(1) 
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at s, 

z{s) = cos 9{t) dt. (3) 

Substituting into ^ we get 



e[d] = ^cEa^{s)el{s)+gpa{s) cos 9 [t) dt^ ds. (4) 

We adopt dimensionless variables in which s, a and e are measured in the units given by 
the scaUng table: 



Variable s a e 
"iMt L ^ ^ 



Here, V is the total volume of the column. V is related to a{s) by 



cL 

V= I a{s)ds. (5) 



The dimensionless versions of ^ and ([5]) are 



e[e]= I !.-a^{s)e^^{s) + Xa{s) [ cos 6 (t) dt} ds, (6) 



V[a] ^ f 
Jo 



^2 

a{s) ds = l (7) 



Where A = is the dimensionless load per unit volume. 

In order to determine the variational BVP for 6 associated with the energy ([6]), it is con- 
venient to change the order of integration in the second term of equation ([6]). Additional 
rearrangements give 

ila^{s)0^,{s) + X cos e{s) [ a{t) dt\ ds. (8) 
1 2 Jo } 

The variational BVP of ([8]) consists of the ODE 

{a^{s)es{s)) +Xsme{s) [ a{t) dt = (9) 
^ Jo 

in < s < 1 and EC 

61(1) = 0, 0^61, = ats = 0. (10) 

Physically, the BC mean that there is a zero angle at the base of the column (clamped) and 
zero torque at the top (free). Equation ([9]) is in fact a local torque balance. An explanation 
of this is given in Appendix IA.3I 
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2.3 Buckling Load 



Recall that A is the dimensionless load per unit volume. For a column tapering given by a{s), 
there is a critical value Ac such that for A < Ac the only solution of the BVP is 9{s) = 0. For 
A > Ac there are non-trivial solutions as well. This is the buckling phenomenon. For such 
A, a small perturbation of the column from the undeflected state can have a smaller energy 
than the undeflected one, therefore the zero solution is no longer an energy minimum. The 
energy minimum is one of the non-zero solutions. For A slightly greater than Ac the energy 
minimizing solution 9{s) will be small. In this case we can study the linearization of ([9]) 
about 6 = 0. The linearization of the BVP ()9|10p is a Sturm-Liouville eigenvalue problem: 

(a^ 9s) + Xc9 [ a{t) dt = 0, in < s < 1 (11) 
^ Jo 

61(1) = a^9s = at s = 0. (12) 

The following analysis assumes, as in [5], that the Sturm-Liouville problem determines a 
discrete sequence of eigenvalues A, and corresponding eigenfunctions which are determined 
up to a multiplicative constant. 



2.4 Eigenvalue Maximization 



What does the tallest column look like? That is, what is the distribution a{s) of cross- 
sectional area which maximizes the smallest eigenvalue Ac of the linearized BVP (jlip subject 
to the volume constraint ([7])? 

Since we have no other A but Ac we will drop the subscript and use A to mean the lowest 
eigenvalue of the linearized BVP. A = X[a] is a functional of a(s). To maximize X[a] subject 
to the volume constraint, we find the variational derivative |^ and solve for the function 
a(s) so that 

Here, /i is a Lagrange multiplier associated with the volume constraint ([7]). 

To find 1^, we introduce a small variation in the cross-sectional area, 6a{s). As a result, A 
and 9{s) will have corresponding small variations SX and 69 (s) respectively. The lineariza- 
tion of equation (fTTI) in 6a, 69 and 6X about a solution a,9,X, results in equation (fT4|) . To 
reduce clutter, the (s) notation is dropped whenever the variable in the parenthesis is s, 
so a, 6a, 9, 69 mean a{s),6a{s),9{s),69{s) respectively. Whenever other variables are used 
they are written explicitly. 

{2a6a9s+a^69s)^ + 6X9 I a{t)dt + X69 a{t)dt + X9 / 6a{t)dt = 0. (14) 

Jo Jo Jo 

The BC for 6a, 69 are 

<56l(l) = (2a6a9s+a'^69s) =0 at s = 0. 
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This is a linear inhomogeneous BVP for 59. Its solvability condition determines the func- 
tional derivative To find the solvability condition, we multiply both sides of (jl4p by 6 
and integrate over < s < 1. Integration by parts, use of BC and rearrangement, leads to 

^ ^{a^es)^ + \e a{t) dt^ 56 ds = !^2a5a9^, - 5X6^ a{t)dt-\5aj^ e'^{t)dt^ds. 



The LHS is zero due to (|lip . and we get an equation relating 5a and 5\: 

S.2ael-x[ e^{t)dt\5ads = 5X [ e^{r) [ a{t)dtdr. (15) 

I Js J Jo Jo 

Here the integration variable in the last term on the right was changed from s to r for added 
clarity. In functional derivative form, equation (jlSj) becomes 

5x _ 2a el -Xj^e^{t)dt 

^~ I,'e^{r)l'a{t)dtdr 

Substituting this into (jl3p results in an integro-differential equation for the maximizing 
cross-sectional area function a and the critical load A: 

e^{t)dt = fi 9^{r) j a{t)dtdr. (16) 
Jo Jo 

Equation (|16p must be solved together with (jlll I12p to produce a(s),0(s), and A. A short 
calculation (see Appendix IA.2p shows that /i = A. 

The integro-differential system (jll|, [T2l [T6|) is still not in a form that we can easily solve. It 
is convenient to convert it into a system of ODE with BC. Differentiating (jl6p with respect 
to s gives the ODE 

2{ael)^ + xe'^ = 

In order to transform equation (jlip into a ODE and get rid of the integral, the variable b{s) 
is introduced: 



h= / a{t)dt. (17) 
Jo 

Differentiating p7p gives a BVP for h{s): 

hs = a, 5(0) = 0. (18) 



Physically, h is the amount of volume above the point s. 
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3 Solution of the BVP 



The system that we want to solve is 



2{ae'l)s + X9^ 



0, 
0, 
0, 



(19) 
(20) 
(21) 



bs - a 



with boundary conditions 



6(1) = 0, 
Os = at s = 0, 



6(1) = 1, 
6(0) = 0. 



(22) 
(23) 



Heuristically, the optimal column is expected to taper to a point as s approaches 0, so a{s) 
is expected to go to zero as s ^ 0. In addition, 6(0) = by the BC (jlSp . Since a and 6 
both multiply equation (jl9p . one expects severe difficulties in the numerical solution as s 
approaches 0, and this is indeed the case. 

This motivates an analysis of the asymptotic structure of the solution to (jl9H23|) as s ap- 
proaches 0. We need to peel away the singularity before attempting any numerical solution. 
In the following sections we use a similarity solution to remove this singularity, transforming 
the equations so that they can be solved numerically. 

3.1 Similarity Solution 

Equations (|19H2ip have a similarity solution. (For information on similarity solutions see, 
for example, [1] or [2].) To find it we study the scaling relations between the variables: 
Specifically, let A, B, S and T be "units" of a, 6, s and 6 respectively. A balance of units in 
ODE'S (fTM 22|) gives the relationships 



T cancels out, and we are left with three equations for A, B and S. In fact, one equation 
is redundant, and we get a simple relation between A, B and S: 



Hence, it is expected that the ODE's are invariant under a scaling transformation, so that if 
9{s),a{s),b{s) are solutions, then so are ^(s/5), 5^a(s/5), 5^6(s/5) for any 5. This suggests 



A^TS-"^ = TB, 
AT^S~^ = r^. 



BS 



A. 



B = S\ 
A = S'^. 
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a similarity solution of the ODE with 

d{s) = aos^ (24) 
b{s) = 5os^. (25) 

We also look for a similar behavior for 9: 

e{s) = BosP. 

Since the equations are linear in 9, we can choose 9q = 1, therefore 

9{s) = sP. (26) 

Balancing units was not enough to determine the exponent of the similarity solution for 9, 
nor for finding the constants aoi ^0- To find these we use ODE's (jl9l — I2ip . 



Substituting (|24l [25]) in equation (pT]) gives the relation 



bo = |. (27) 
Further substitution of (piH271) into (|19|20p yields 

p{p + 5) = -7, 
^/(2p+i) = _^, 

where 7 = Setting the two expressions for —7 equal gives a polynomial equation for p: 

2p{p + 5) =p2(2p + l). 

The solutions to this equation are 0, — 2, |. Each solution for p implies a solution for 7. By 
looking at equation (jl6p we see that 7 must be non- negative. Indeed, when s = 1, the LHS 
is non-negative, and /i = A implies A > 0. Since oq is positive, 7 > 0. The values 0, —2, | 
for p imply the values 0, 6, —75 for 7, respectively. Therefore, the only admissible values for 
p are 0, —2. p = implies 7 = 0, and this is physically non-interesting: When 7 = 0, A = 0, 
hence the column has no weight. Either there is no gravity, or the density of the material 
is zero. We want to study the first positive eigenvalue A > 0. The solution p = —2 gives 
7 = 6 hence, oq = The similarity solution is found to be 

a{s) = (28) 

ks) = (29) 
9{s) = (30) 
It is easy to see that the similarity solution (j28H30p satisfies the BC (j23p as s ^ 0. 
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3.2 Peeling Away the Singularity 



While the similarity solution (j28H30p has the correct asymptotic behavior as s ^ 0, it does 
not satisfy BC (f22l) at the base s = 1 (because 9{l) ^ 0). 

We use the similarity solution to find a solution of the ODE that satisfies the BC at both 
ends. Any solution can be written as a product of the similarity solutions and new dependent 
variables a, /3, r: 

a = aa, (31) 
b = b(3, (32) 
e = eT. (33) 

Substituting these expressions for a, b and 6 in the differential equations ()19H2ip results in 
ODE's for a, (3 and r. These ODE's are homogeneous in s, hence we use t = — Ins as the 
independent variable. The ODE for a(t), /3(t), r(t) are the autonomous system (AS): 

(3-L>) (a2(2 + D)r) -6/3r = 0, (34) 
(3 + Z?) (a [(2 + Z))r]2) - 12r2 = 0, (35) 
(4 - L>) /3 - 4a = 0. (36) 

In (j34H36p . D denotes differentiation with respect to the variable t. The BC ()22l23p are 
also transformed. At the base, s = 1, t = the base BC ()22p transform into 

m = y , r(0) = 0, (37) 



At the tip, as s ^ 0, t ^ oo, the BC (j23j) transform into 

r(i) ^1 as i — > oo, f3{t) -^1 as t ^ cx). (38) 

The ODE's above, although written in an implicit form, can be rewritten as explicit ex- 
pressions for Ttt,fit,ctt as functions of T,Tt,(3,a. Note that due to the specific choice of 
dependent variables, A is no longer a parameter in the AS. It only appears in the BC at the 
base. It is easy to see that a = /3 = l,r = lisa critical point of the AS (pH - [36|) . Setting 
a = /3 = T = lin ([33H32|) recovers the similarity solution (p8H30|) . 

We seek a solution that satisfies the BC at t = and approaches the critical point as t ^ oo. 
This solution satisfy all the BC. A solution that approaches a critical point as its limit at 
t = oo is a solution on the stable manifold of the critical point. 

Since the critical point has both a stable an unstable manifold, it is difficult to start on 
the stable manifold far from the critical point and numerically follow the solution into the 
critical point. The smallest numerical error will cause the solution to deviate on a growing 
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solution that does not approach the critical point. Instead, we identify the stable manifold, 
and follow the solution backwards in t from the critical point until the BC at t = are 
satisfied. This way, numerical errors will not be amplified. In fact, numerical errors will 
either decay or remain constant when solving in this direction. The following section shows 
that the stable manifold is one dimensional. Thus, there is only one direction in which to 
start the solution near the critical point. Solving the AS backwards from the critical point 
will follow the stable manifold. Since the BC at the base define a plane of co-dimension 1, 
the stable manifold is expected to intersect this plane. Once the intersection point is found, 
the solution is fully determined. 



3.3 Stable Manifold 



To identify the stable and unstable manifolds of the critical point, we linearize the AS 
around the critical point(l, 1, 1), resulting in the linear ODE 



/D{D + 5) 
\0~ 



D + 3 

6 4{D - 3) 
(L> - 4) 4 




(39) 



Here, {6T,5P,6a) are the deviations from (1,1,1). (j39p is a linear system of ODE's, we 
therefore look for solutions of the form 

{6T,6/3,6a) = {6To,5Po,5ao)e'^^ . 



Substituting this form of solution into equation ()39l) yields a generalized linear eigenvalue 
problem: 



^q{q + b) q + 3 

q{q - 1) 6 A{q - 3) 
(g - 4) 4 




(40) 



Equation (j40p has non-trivial solutions for 4 different values of g : gi = 0, ^2 = Ij^Zs ~ —5.5208, 
and g4 ~ 6.5208. The solutions are spanned by the vectors given in the table below. As we 





Si 


^2 


S3 


Si 




1 


-2 


0.876733 


-0.126733 







3 


0.420133 


-1.5868 







4 


1 


1 



Table 1: possible solutions to equation ([iOl) 

are looking for the stable manifold, the only solution of interest is the one with a negative q, 
i.e. solution (93,53) Readers noticing the conspicuous form of solution (92, 5*2) are referred 
to Appendix lA.ll for brief discussion of it. 
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To find the numerical solution of (j34H36p . we start near the stationary point, on the line 
tangent to the stable manifold, and solve the AS backwards in t. The BC (|37p at the base 
define a surface B in the r, Tt, (3,a phase space: 



The first condition, r = 0, is used as a stopping condition the base of the column. The 
second, /? = identifies A once the base is found. Once A is found, the original functions, 
9,a,b are known from ([HTH55|) . 

3.4 Numerical results 

Using Matlab 's ODE solver ode45, which uses the fourth-order Runge-Kutta method, the 
AS was solved. The solver was used with relative error tolerance of 10~^ and an absolute 
error tolerance of 10~^. The initial value for (r, rj,a,/3) was 



The values for To,aO)/5o and q are taken from ((j'3,5'3). Since the stable manifold is one 
dimensional there are two distinct options for starting near the critical point. In terms of 
(j4ip the two options are choosing 6 to be positive or negative. The solution with positive 
6 fails to satisfy the BC at the base. Geometrically, the trajectory is going in the wrong 
direction in the phase space and does not intersect the plane B. With 6 = —0.0001 the BC 
at the base are satisfied (that is, r vanishes) at At = —1.7114. See figure [3l 

Since at the base, A = the value for A is given from the value of /5. By this formula 
A = 134.1944. Higher accuracy can be achieved by starting closer to the critical point 
(smaller 5), and setting the solver to a smaller tolerance. 

Equations (f3TH33l) uniquely determine 9, a and b given a, (3, r and A. Figure [3] shows the 
resulting cross-section a(s) of the A-maximizing column tapering. 

4 Other Boundary Conditions 

The fact that A appears only in the BC at the base has an interesting implication for this 
problem: If we want to solve the problem with a different BC at the base, we will have 
the same analysis at the tip and the same similarity solution. Since the stable manifold is 
one dimensional, we will have to start in the same way while solving backwards. The only 
difference will be when to stop. 

For instance, consider the problem of the tallest hinged column. A hinged column is a 
column that instead of being constrained to be vertical at the base, is forced to have zero 
torque at the base. 



B = \ (T,Tt,i5,a) r = and — - /? = 




{T,Tt,a,(3) = (1,0,1,1) + 6 ■ (ro,gTo,ao,/3o)- 



(41) 
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0.5 1 1.5 0.5 1 

t S 

(a) (b) 

Figure 3: The solutions of the tallest clamped column. Figure [3)^a) shows T{t). It satisfies 
r(0) = and as t increases it approaches 1. Figure [3l|b) shows a{s) (solid line) and the 
similarity solution (dashed line). Since our solution for a{t) extends only up to a finite t, it 
isn't drawn for small values of s. For small s (i.e. large t) we can use the similarity solution 
to extend the solution to s = 0. 



Note: The lowest eigenvalue of the hinged column is zero. This is because at no gravity 
the column can be set at any angle. To find the first non-trivial solution it would appear 
that we need to add an additional orthogonality condition to make sure that we get the 
next eigenfunction. A short calculation shows that the boundary conditions and the ODE 
already guarantee the orthogonality of solutions to the trivial, constant-angle, solution. 

Formally, the difference is in the BC. In (llOp . ^(1) = is replaced by the zero torque 
condition 

a'^0s = O at s = l. (42) 

The translation of this BC to r, a, (3 is 

(2t + Tt) = att = 0. 

The solution to this problem would be the design of the tallest column that is hinged at the 
base and free at the top. Using the same numerical method and the same initial conditions 
as before, a numerical solution was found. After At = —1.9470, the BC, a'^{2T + Tt) = 0, is 
satisfied. 

Again, since A = we have the value of A. For these BC, the value is A = 222.7366. The 
cross-section and "peeled" torque are shown in Figured! 

The results predict a larger eigenvalue for the building with hinged base than with a clamped 
base. This might be non-intuitive until one notices that the hinged column has within it a 
shorter clamped column. This column is also buckling and since it is shorter, its buckling 
load will be larger than the full length clamped beam. This is also the case with simple 
loaded beams. The hinged loaded beam will be able to support a larger load than the 
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Figure 4: The solution of the tahest hinged column. Figure IHa) is the "peeled" torque, 
a^(2r + rj), as a function of t. It goes to zero at t = and since a, r go to the critical point 
as t increases, it goes to 2. Figure IHb) shows a{s) (solid line) and the similarity solution 
(dashed line). As before, a(s) isn't drawn for small values of s. For small s (i.e. large t) we 
can use the similarity solution to extend the solution to s = 0. 

clamped one. The reason for this non-intuitive result is the assumption that the column 
actually remains vertical. 

5 Discussion 

The derivation of the equations above relies on the use of a variational derivative in order 
to find the cross-section a(s) that produces a stationary eigenvalue A. For this procedure 
to be valid, the operator, BVP ([H [TO]) should have discrete eigenvalues. As Cox and 
McCarthy have already shown, this is not always true. Nevertheless, the Keller and Niordson 
solution is very plausible. As shown in [5], adding a small weight to the top of the column 
and letting the weight tend to zero produces a sequence of optimization problems. These 
problems converge to the self-weight problem, and the solutions (which in this case must 
exist) converge to the described solution. 

We confined ourselves to demonstrating the use of the self similar solution in order to 
simplify the numerical solution of the given equations. Regardless of their derivation, the 
equations themselves require a solution and the method of peeling away the singularity 
using the similarity solution does away with the need of a iteration scheme. This allows for 
any standard ODE solver to tackle the remaining problem. This method can be applied to 
other problems that have singularities. In a companion paper we employ a similar analysis 
to find the tapering of a javelin so that its lowest vibration mode has the highest possible 
frequency. 
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A Appendix 



A.l Specific Solutions in the Linearized System 

We turn our attention to the conspicuous value of the eigenvalue ^2 = 1 and respective 
eigenvector S2 = (3,4,-2) which matches exactly with the exponents of the similarity 
solution. Here we show that this solution of the linearized ODE is to be expected. 



Let 



where 



i = f (x) 



(43) 



X 



be an ODE. We further assume that 



X(s) 



\ Xn J 



a2S 



P2 



aj / 



(44) 



is a solution to the ODE. Any other solution Y can be written as 

Y{s) = A{s)y.{s) 



(45) 



where A{s) is a diagonal matrix. The ODE that A{s) must satisfy can be found by substi- 
tuting (I45|44p into (03]). 

f(A(s)X(s)) = f (Y(s)) = Y(s) = i(s)X(s) + A{s)±{s) = i(s)X(s) + A{s)Py.{s)s-^ . 

Here, P is a diagonal matrix with Pa = pi. In other words, A{s) must satisfy 

f(A(s)X(s)) = i(s)X(s) + A{s)PX{s)s-\ (46) 

Note that A{s) = / is a solution of this system. Therefore, / is a critical point of the 
system for A. The stable and unstable manifolds of this system are important for the type 
of solution method described in this paper. The following proposition shows the existence 
of a specific unstable direction about the critical point /. 



Proposition. The linearization of (|46p around I has a solution 
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In the scaled variable, t = 
P. 



6A{s) = Ps-\ (47) 
— Ins, this corresponds to an eigenvalue of 1 and an eigenvector 



Proof. The linearization of ()46p is 

JF 6A{s)X{s) = 6'A{s)X{s) + 6A{s)PX{s)s~\ (48) 

Here, dA is a diagonal matrix (the perturbation from the solution I) and JF is the Jacobian 
matrix of f evaluated at ^(s)X(s). To find out more about JF, we differentiate the equality 

Ps-^X = f (X) (49) 

with respect to s: 

P{P - I)X{s)s-'^ = J£P:S.{s)s-^. (50) 

Equation (j49p is derived by substituting (j44p into (|43p . To show that 6A{s) = Ps~^ solves 
equation (j48]), we substitute (|47|) into (f48l) : 



-Ps-'^X{s) + P^X{s)s"'^ = JfPX{s)s-\ 



From equation (|50l) we substitute for the RHS 

-Ps-2x(s) + p2x(s)s-2 = p(p - /)X(s)s~2 
This shows that indeed 6A{s) = Ps~^ is a solution of the linearized ODE. □ 

A. 2 Calculation of the Lagrange multiplier 

To calculate the value of /x in (fT6p . we multiply the equation by a(s) and integrate from 
to 1: 

1 /•! pi rl rr rl 

2a^elds-\j a I 6"^ (t) dt ds = fi 0^{r) a{t)dtdr ads. 
^0 Js Jo Jo Jo 

Since the RHS of (fTUp is a constant and a{s) ds = 1, the RHS stays unchanged. In the 
LHS, integrating the first term by parts, changing the order of integration of the second 
term and using the BC gives 



-2 / {a^9s)jds - \ [ 9^ [ 
Jo Jo Jo 



1 fs pi pr 

a{t)dtds=fjL I e'^{r) / a(t)dtdr. 
Jo Jo 
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Lastly, using ()TT|) we can see that 



f'l rs f'l rs f'l rr 

2A / a{t)dtds-\j a{t)dtds = fi e^{r) a{t)dtdr 

Jo Jo Jo Jo Jo Jo 

Therefore, A = /i. 



A. 3 Torque Balance 




Figure 5: The torque at point s is due to the sum of the torques from points s' < s along 
the column. Each point creates a torque a(s') \Axds' (recall that A is the non-dimensional 
gravity). The curvature at the point s is proportional to this torque. 



The differential equation ([9]) is a local torque balance. Integrating equation ([9]) and using 
the second of BC gives 



(51) 



a'^{s)es{s) = X sme{s') / a{t)dtds. 
Jo Jo 

Changing the order of integration in the RHS of (I5ip we see that 

a{s') I sin^(t) dtds'. 



(52) 



The inner integral is the horizontal displacement Ax of the point s' from the point s. The 
sections of the column above s exert a gravity generated torque on the point s. This is the 
RHS of (f52]) . The LHS expresses the elastic response to this torque: The curvature Og of 
the beam at s is proportional to the imposed torque. See figure ([5|). 



17 



References 

[1] G. I. Barenblatt, Scaling, self- similarity, and intermediate asymptotics, Cambridge Uni- 
versity Press, New York, 1996. 

[2] G. W. Bluman and S. C. Anco, Symmetry and integration methods for differential equa- 
tions, Applied mathematical sciences, vol. 154, Springer- Verlag, 2002. 

[3] Steven J. Cox and C Maeve McCarthy, The shape of the tallest column, SIAM Jom'nal 
of Apphed Math Anals 29 (1998), no. 3, 547-554. 

[4] J. B. Keller, The shape of the strongest column. Arch. Rational Mech. Anal 5 (1960), 
275-285. 

[5] J. B. Keller and F. I. Niordson, The tallest column, J. Math. Mech. 16 (1966), 433-446. 

[6] Philip G. Kirmser and Kuo-Kuang Hu, The shape of the ideal column reconsidered. The 
Mathematical Intelligencer 15 (1993), no. 3, 62-68. 

[7] C. Maeve McCarthy, The tallest column — optimality revisited. Journal of computa- 
tional and applied mathematics (1999), no. 101, 27-37. 

[8] Charles A. Stuart, Buckling of a heavy tapered rod. Journal de Mathmatiques Pures et 
Appliques 80 (2001), no. 3, 281-337. 



18 



